Numerical Simulations of the Low-Velocity Impact Response of Semicylindrical Woven Composite Shells

This paper presents an efficient and reliable approach to study the low-velocity impact response of woven composite shells using 3D finite element models that account for the physical intralaminar and interlaminar progressive damage. The authors’ previous work on the experimental assessment of the effect of thickness on the impact response of semicylindrical composite laminated shells served as the basis for this paper. Therefore, the finite element models were put to the test in comparison to the experimental findings. A good agreement was obtained between the numerical predictions and experimental data for the load and energy histories as well as for the maximum impact load, maximum displacement, and contact time. The use of the mass-scaling technique was successfully implemented, reducing considerably the computing cost of the solutions. The maximum load, maximum displacement, and contact time are negligibly affected by the choice of finite element mesh discretization. However, it has an impact on the initiation and progression of interlaminar damage. Therefore, to accurately compute delamination, its correct definition is of upmost importance. The validation of these finite element models opens the possibility for further numerical studies on of woven composite shells and enables shortening the time and expenses associated with the experimental testing.


Introduction
Due to its distinctive combination of high strength, low weight, and exceptional fatigue resistance, composite materials have grown in popularity. Nevertheless, lowvelocity impacts that may happen during handling, transit, maintenance, and service might harm them. These collisions may result in localized structural damage, which may diminish the composite material's strength, stiffness, and durability. To ensure the dependability and safety of composite structures, it is crucial to comprehend the behavior of composite materials under low-velocity impacts.
In a low-velocity impact event, composite materials undergo various stages of damage. Firstly, when the impactor contacts with the material, a sudden increase in stress and strains is experienced in a very localized region. Secondly, microcracks start to develop in the material matrix, which then propagate through the material and spread to the adjacent laminas and neighboring interface regions. At this stage, the damage progression causes the separation of the adjacent laminas of the composite material, which is known as delamination. This damage mechanism, alongside with matrix cracking, fiber failure and perforation, can significantly reduce the strength and stiffness of the material. The extent and severity of damage depends on the properties of the composite material, the impactor's properties, and the impact energy. The development of reliable models capable of predicting these damage stages is of highest importance for the analysis and design of composite structures that are resistant to low-velocity impacts.
In this context, the behavior of composite structures has been extensively studied using finite element (FE) models in a variety of industries, including the marine, automotive, and aerospace sectors [1][2][3][4][5]. The impact response of composite flat plates has been extensively studied using numerical models to examine the impacts of various parameters such as material characteristics, impact energy, and geometry. Yet, there are relatively few numerical studies that analyze the impact dynamics on cylindrical shells, particularly when it comes to composites made of woven fabric. The most relevant numerical studies that were conducted on cylindrical shells include the work developed by Kim et al. [6] in which they observed that the contact force increases with the curvature of shell-shaped composite laminates, while the deflection and contact time decrease. It was also observed that the impactor's velocity has a greater influence on the contact force than the impactor's mass, which is similar to the impact response of composite plates [6][7][8]. This is justified by the fact that the kinetic energy of the impactor increases linearly with the increase in mass and quadratically with velocity. The contact force and deflection histories for composite laminated cylindrical shells with convex and concave shapes were analyzed by Choi [9]. The author found the same contact force and central deflection histories for both shapes. Kistler and Wass [10] performed a numerical study on unidirectional (UD) laminated cylindrical shells and identified scaling relationships between impact energy, momentum, mass, and velocity, while Zhao and Cho [11] investigated the impact-induced damage initiation and propagation of UD composite shells and found that the damage propagates differently from composite flat plates. Another study was performed by Kumar et al. [12] to study the impact response and impact-induced damages of cylindrical UD composite laminate shells using a 3D finite element formulation. Recently, Khalili et al. [13] used FE analysis to investigate the impact response of UD composite laminate plates and shells structures under low-velocity impact loads and optimize the procedure for future work. Albayrak et al. [14] conducted an experimental and numerical investigation to study the geometrical effect on low-velocity impact behavior for curved composites with a rubber interlayer. They found that the curved surface geometry affects the absorbed energy and that increasing the width of the laminate while keeping the height constant results in higher impact energy absorption.
Overall, these studies provide valuable insights into the behavior of composite laminate structures under low-velocity impacts and can inform about the design and optimization of such structures for various applications. Nevertheless, none of these numerical studies was dedicated to the development of FE models capable of analyzing the lowvelocity impact response of semicylindrical woven fabric composite shells. Therefore, the main goal of this paper is to develop reliable and efficient FE models capable of predicting the impact behavior of these materials. For this purpose, constitutive modes that consider the intralaminar and interlaminar progressive damage were implemented in the explicit finite element approach using ABAQUS/Explicit [15]. The validation of the FE models was carried out using the authors' previous work on the experimental assessment of the effect of thickness on the multi-impact response of semicylindrical composite laminated shells [16]. To facilitate the understanding, the nomenclature used in this paper is listed in the Nomenclature.

Material and Experimental Procedure
Composite semicylindrical shells were produced using a matrix based on an AROPOL FS 1962 polyester resin and a MEKP-50 hardener (both supplied by SF Composites, Mauguio, France). A bi-directional E-glass woven fabric (taffeta with 210 g/m 2 ) was used as reinforcement, and the composite was produced by hand lay-up with 9 woven fabric layers (corresponding to 1.6 mm of final thickness). In order to ensure a constant fiber volume fraction and uniform thickness, as well as to eliminate any air bubbles, the laminates were placed inside a vacuum bag immediately after impregnation. The manufacturing process culminated with curing at 40 • C for 24 h. More details about the materials and manufacturing process can be found in [16][17][18]. Figure 1 shows the specimens' dimensions and the schematic view of the test conditions. volume fraction and uniform thickness, as well as to eliminate any air bubbles, the laminates were placed inside a vacuum bag immediately after impregnation. The manufacturing process culminated with curing at 40 °C for 24 h. More details about the materials and manufacturing process can be found in [16][17][18]. Figure 1 shows the specimens' dimensions and the schematic view of the test conditions. Finally, low-velocity impact tests were performed on a drop weight testing machine IMATEK-IM10 (Old Knebworth, UK), which is described in detail in [19]. These tests were carried out according to ASTM D7136 standard, at room temperature, and using an impactor diameter of 10 mm with a mass of 2.826 kg. The energy of 5 J was used to promote visible damage but without full perforation. More details about these tests can be found in [16].

Damage Models
Two constitutive models were used in this study to simulate the material damage caused by low-velocity impact loads in semicylindrical composite laminate shells: (i) a continuum damage model (CDM) at the lamina level to account for intralaminar damage; and (ii) a surface-based cohesive damage model (S-BCM) at the lamina interface to account for interlaminar damage.
The built-in constitutive model for fabric-reinforced composites available in Abaqus/Explicit [15], developed by Johnson [20] and based on Ladeveze and Ledantec work [21], was used to evaluate the complex damage progression at the intralaminar level. When the user-defined material is named with the string "ABQ PLY FABRIC," this model can be used as a built-in VUMAT user subroutine [22].
The maximum stress criterion determines the damage initiation, and the fracture energies serve as the basis for the damage evolution model, which regulates the decline in stiffness. In this way, the following damage activation functions, and , are used to compute the elastic domain at any given time, Finally, low-velocity impact tests were performed on a drop weight testing machine IMATEK-IM10 (Old Knebworth, UK), which is described in detail in [19]. These tests were carried out according to ASTM D7136 standard, at room temperature, and using an impactor diameter of 10 mm with a mass of 2.826 kg. The energy of 5 J was used to promote visible damage but without full perforation. More details about these tests can be found in [16].

Damage Models
Two constitutive models were used in this study to simulate the material damage caused by low-velocity impact loads in semicylindrical composite laminate shells: (i) a continuum damage model (CDM) at the lamina level to account for intralaminar damage; and (ii) a surface-based cohesive damage model (S-BCM) at the lamina interface to account for interlaminar damage.
The built-in constitutive model for fabric-reinforced composites available in Abaqus/Explicit [15], developed by Johnson [20] and based on Ladeveze and Ledantec work [21], was used to evaluate the complex damage progression at the intralaminar level. When the user-defined material is named with the string "ABQ PLY FABRIC", this model can be used as a built-in VUMAT user subroutine [22].
The maximum stress criterion determines the damage initiation, and the fracture energies serve as the basis for the damage evolution model, which regulates the decline in stiffness. In this way, the following damage activation functions, F α and F 12 , are used to compute the elastic domain at any given time, where σ α and σ 12 are the effective normal and shear stresses, respectively, X α is the tensile/compressive strength, S 12 is the shear strength, and r α and r 12 are the corresponding damage thresholds, which are initially set to 1. Once damage is predicted, the elastic-stressstrain relations are given by, where d 1 and d 2 are the damage variables associated with fiber fracture along directions 1 and 2, respectively, and d 12 is the damage variable associated with the matrix microcracking due to shear deformation. These variables are determined with Equations (4) and (5), where r α and r 12 stand for the damage thresholds for axial and shear loads, respectively, G α f stands for the fracture energies, g α 0 to the elastic energy density, L e stands for the characteristic length of the element, and α 12 stands for the shear damage parameter.
The non-linear behavior of the matrix, due to microcracking, dominates the shear damage response at the intralaminar level and includes both stiffness reduction and plasticity. The latter is defined with the following yield and hardening functions; see Equations (6) and (7), respectively.
where σ y0 and σ 0 are the initial effective shear stress and shear yield stress, ε pl is the plastic strain due to shear deformation, and C and the superscript p correspond to the coefficient and power term in the hardening function, respectively. Using the two-step homogenization methodology described by Liu et al. [23], the stiffness properties of the woven fabric composite laminas were estimated from the constituents' properties of the tested specimens in order to validate the numerical model based on the experimental evidence presented in [16]. The values of the remaining intralaminar material properties were taken from the literature. In this way, the fracture toughness values, the damage evolution parameters, α 12 and d max 12 , and the shear plasticity parameters, σ y0 , C and p, were taken from impact studies on E-glass laminates [24][25][26], whereas the strength properties were taken from impact studies on woven E-glass/polyester composite laminates [27][28][29]. The intralaminar material parameters needed to specify the material model in the VUMAT subroutine are shown in Table 1. It is noteworthy to mention that a preliminary parametric study was performed using a coarser FE mesh to assess how a reasonable variation of the shear plasticity parameters could impact the results. It was found that they have a negligible effect on the numerical solutions. Regarding the strength properties, the data found in the literature vary substantially, depending on the manufacturing process, type e-glass fiber, volume fraction, etc. Taking this fact into consideration, the values employed are based on averaged values and again, a preliminary parametric study was performed with a coarser FE mesh to assess which values best fit the experimental evidence. Finally, to account for the interlaminar damage, the bond between the laminas of the composite laminate was modeled using cohesive surfaces. This approach is primarily intended for negligible small interface thicknesses and offers very similar capabilities to cohesive elements. The cohesive behavior is defined as a surface interaction property, and identically to the cohesive elements, it is governed by a traction-separation constitutive model. The properties employed in the surfaced-based cohesive model shown in Table 2 were extracted from the literature [30][31][32][33][34][35][36][37][38]. Table 1. Intralaminar properties.

Property Symbol Units Value
Stiffness properties k n = k s = k t N/mm 3 10 6 Strength properties τ 0

MPa 30
Fracture toughness The cohesive stiffness in this study is set at 10 6 N/mm 3 , as suggested by Camanho et al. [30]. In addition, it is considered that its value is the same for all directions, that is, k n = k s = k t , as used in [31][32][33] with satisfactory results. It is noteworthy to mention that considering high values for the cohesive stiffness potentially results in convergence problems. On the other hand, the use of low values may affect the global stiffness and thus compromise the validation of the FE model [32]. A value of η = 1.45 was considered for the interaction parameter in the definition of the cohesive model [34][35][36][37][38]. Identically to the intralaminar properties, there is a wide range in the data for the interlaminar strength parameters and fracture toughness. Given this information, the values used are averages, and preliminary parametric analysis using a coarser FE mesh was completed to determine which values the best-matched experimental data.

Finite Element Model
The FE model was created using the ABAQUS/Explicit FE code [15] taking into account the dimensions of the nine-layer laminate specimens evaluated in [16]. The tested specimens had a semicircular internal radius of 50 mm, a length of 100 mm and an average thickness of 1.6 mm ( Figure 2).  To replicate the experiments, two fixed rigid body supports (a lateral and a bottom support) were included in the FE model. Only one-fourth of the semicylindrical composite laminate was generated, taking use of the geometric symmetry of the model to reduce the computing cost of the numerical simulations. The yz-plane face = = = 0 and one of the xy-plane faces = = = 0 were therefore added to the symmetry boundary conditions. The impactor was modeled with a lumped mass fixed on a reference point at its center of mass equivalent to the experiments and with a hemispherical head with a diameter of 10 mm. Only the displacements in the y-direction were permitted = = 0 , and all the rotations of the impactor were constrained , , = 0 . Each lamina was discretized with SC8R continuum shell elements (eight-node hexahedron) with reduced integration and stiffness hourglass formulation. The orientations of the materials along the semicircular cross-section were taken into account when defining the local coordinate system of the laminas. R3D4 discrete rigid elements were used to model the impactor. The lamina was modeled with cohesive surfaces; thus, no element specification was required.
The surface-to-surface contacts between the composite laminate, the metal impactor, and the metal supports were simulated using the penalty enforcement contact method from Abaqus/Explicit [15]. In the interface of the composite laminas, which experience friction after being entirely delaminated, this contact formulation was also defined. The friction coefficient values, , used in this work for fully damaged interfaces and metalcomposite contacts were taken from [39,40]. A value of = 0.3 was specified for the contact between the metal hemispherical head of the impactor and the upper surface of the composite laminate, and a value of = 0.7 was specified for the contact between the metal surfaces of the supports and the composite laminate surfaces. For the interfaces, a value of = 0.5 was considered.

Numerical Results
Several numerical simulations were performed to determine the influence of the FE mesh discretization and mass scaling on the efficiency and reliability of the FE model. To be able to analyze how these parameters affect the numerical predictions, the most relevant load and energy histories are presented, as well as the maximum force, displacement, and contact time.
The use of cohesive surfaces implies that its elements' characteristic length matches the characteristic length of the continuum shell elements defined for the laminas. Given that the interlaminar damage surface-based cohesive damage model implementation To replicate the experiments, two fixed rigid body supports (a lateral and a bottom support) were included in the FE model. Only one-fourth of the semicylindrical composite laminate was generated, taking use of the geometric symmetry of the model to reduce the computing cost of the numerical simulations. The yz-plane face U x = R y = R z = 0 and one of the xy-plane faces U z = R x = R y = 0 were therefore added to the symmetry boundary conditions. The impactor was modeled with a lumped mass fixed on a reference point at its center of mass equivalent to the experiments and with a hemispherical head with a diameter of 10 mm. Only the displacements in the y-direction were permitted (U x = U z = 0), and all the rotations of the impactor were constrained R x,y,z = 0 .
Each lamina was discretized with SC8R continuum shell elements (eight-node hexahedron) with reduced integration and stiffness hourglass formulation. The orientations of the materials along the semicircular cross-section were taken into account when defining the local coordinate system of the laminas. R3D4 discrete rigid elements were used to model the impactor. The lamina was modeled with cohesive surfaces; thus, no element specification was required.
The surface-to-surface contacts between the composite laminate, the metal impactor, and the metal supports were simulated using the penalty enforcement contact method from Abaqus/Explicit [15]. In the interface of the composite laminas, which experience friction after being entirely delaminated, this contact formulation was also defined. The friction coefficient values, µ, used in this work for fully damaged interfaces and metal-composite contacts were taken from [39,40]. A value of µ = 0.3 was specified for the contact between the metal hemispherical head of the impactor and the upper surface of the composite laminate, and a value of µ = 0.7 was specified for the contact between the metal surfaces of the supports and the composite laminate surfaces. For the interfaces, a value of µ = 0.5 was considered.

Numerical Results
Several numerical simulations were performed to determine the influence of the FE mesh discretization and mass scaling on the efficiency and reliability of the FE model. To be able to analyze how these parameters affect the numerical predictions, the most relevant load and energy histories are presented, as well as the maximum force, displacement, and contact time.
The use of cohesive surfaces implies that its elements' characteristic length matches the characteristic length of the continuum shell elements defined for the laminas. Given that the interlaminar damage surface-based cohesive damage model implementation yields mesh-dependent results, the FE mesh discretization of the laminas needs to be defined based on the characteristic length of the cohesive surface elements. In other words, the FE mesh size of the interface defines the size of the FE mesh employed in the whole model.
To find a balance between the computational cost of the solution and the accurate computation of the fracture toughness in the interlaminar damage model, which results in delamination, a parametric study was conducted to optimize the characteristic length of the elements of the FE mesh. Based on the work of Hilleborg et al. [41], Turon et al. [31] purposed using Equations (8) and (9) where M is a parameter that depends on the cohesive zone model, and N e is the number of elements in the cohesive zone. The lowest value derived from the equations is l e,I I = 0.31 mm, with the assumptions that M = 1, as suggested in [23,33,35], N e = 5 to properly establish the cohesive zone [42][43][44], and the baseline properties of the laminas and of the cohesive zone. Consequently, the baseline FE mesh was generated with l e = 0.3 mm, which corresponds to an aspect ratio of 1.6. Notice that this value is in good agreement with that used by Lopes et al. in [35]. This baseline FE model contains approximately half a million elements and one million nodes. Therefore, to shorten the computing time for the solutions, a semi-automatic mass scaling was uniformly applied to the entire model with a target time increment of 1 × 10 −7 . Notice that mass scaling artificially increases the mass of the structure to reduce the frequency of the dynamic response and allows the time step of the simulation to be increased. In this study, a mass increase of 1.8% was obtained for the baseline FE model, resulting in a reduction in the computational cost of the solutions of about 51%. To assess the impact of mass scaling on the load and energy history curves, the results obtained with the baseline FE model are shown in Figure 3

model.
To find a balance between the computational cost of the solution and the accu computation of the fracture toughness in the interlaminar damage model, which res in delamination, a parametric study was conducted to optimize the characteristic len of the elements of the FE mesh. Based on the work of Hilleborg et al. [41], Turon et al. purposed using Equations (8) and (9) to calculate the characteristic length of the elem in the direction of the crack propagation for fracture modes I, II, and III in orthotro composite materials, is a parameter that depends on the cohesive zone model, and is the num of elements in the cohesive zone. The lowest value derived from the equations is 0.31 mm, with the assumptions that = 1 , as suggested in [23,33,35], = 5 properly establish the cohesive zone [42][43][44], and the baseline properties of the lam and of the cohesive zone. Consequently, the baseline FE mesh was generated with 0.3 mm, which corresponds to an aspect ratio of 1.6. Notice that this value is in good ag ment with that used by Lopes et al. in [35]. This baseline FE model contains approximately half a million elements and one lion nodes. Therefore, to shorten the computing time for the solutions, a semi-autom mass scaling was uniformly applied to the entire model with a target time incremen 1 10 . Notice that mass scaling artificially increases the mass of the structure to red the frequency of the dynamic response and allows the time step of the simulation to increased. In this study, a mass increase of 1.8% was obtained for the baseline FE mo resulting in a reduction in the computational cost of the solutions of about 51%. To as the impact of mass scaling on the load and energy history curves, the results obtai with the baseline FE model are shown in Figure 3     It is possible to observe that these curves include oscillations brought on by the e wave and produced by the models' vibrations [45,46]. The numerical predictions f numerical maximum load, maximum displacement, and contact time are compared ble 3.    It is possible to observe that these curves include oscillations brought on by the wave and produced by the models' vibrations [45,46]. The numerical predictions numerical maximum load, maximum displacement, and contact time are compare ble 3.  It is possible to observe that these curves include oscillations brought on by the elastic wave and produced by the models' vibrations [45,46]. The numerical predictions for the numerical maximum load, maximum displacement, and contact time are compared in Table 3. It can be observed that the impact response of the curves is similar with and without mass scaling. Its application has a negligible effect on the contact time and no impact on the maximum displacement. Only for the maximum force do the numerical predictions differ by 8.8%; its value is higher when mass scaling is employed. Overall, the results indicate that the defined semi-automatic mass-scaling parameters have an acceptable impact on the numerical predictions and significantly lower the computational cost of the solutions.
To determine if a coarser FE mesh discretization would produce a good trade-off between the accuracy of the solution and computing cost, a parametric study was performed using increasingly coarser FE meshes, ranging from l e = 0.3 mm to l e = 2 mm. The FE models employed to analyze the effect of the FE mesh discretization are shown in Figure 6.

FOR PEER REVIEW 9 of 17
It can be observed that the impact response of the curves is similar with and without mass scaling. Its application has a negligible effect on the contact time and no impact on the maximum displacement. Only for the maximum force do the numerical predictions differ by 8.8%; its value is higher when mass scaling is employed. Overall, the results indicate that the defined semi-automatic mass-scaling parameters have an acceptable impact on the numerical predictions and significantly lower the computational cost of the solutions.
To determine if a coarser FE mesh discretization would produce a good trade-off between the accuracy of the solution and computing cost, a parametric study was performed using increasingly coarser FE meshes, ranging from = 0.3 mm to = 2 mm. The FE models employed to analyze the effect of the FE mesh discretization are shown in Figure 6. The numerically predicted force-time and energy-time results are shown in Figure 7 and Figure 8, respectively, and the force-displacement results are shown in Figure 9. The results show that the use of coarser FE meshes ( = 1 mm and = 2 mm) induces higher oscillations on the force-time and force-displacement curves. Nonetheless, the maximum force, maximum displacement, and contact time values are barely affected. This can be observed in Table 4, where the values and percentage difference between = 0.3 mm and the remaining element lengths are presented. The numerically predicted force-time and energy-time results are shown in Figures 7  and 8, respectively, and the force-displacement results are shown in Figure 9. The results show that the use of coarser FE meshes (l e = 1 mm and l e = 2 mm) induces higher oscillations on the force-time and force-displacement curves. Nonetheless, the maximum force, maximum displacement, and contact time values are barely affected. This can be observed in Table 4, where the values and percentage difference between l e = 0.3 mm and the remaining element lengths are presented.
Therefore, it is clear that the FE mesh discretization, within the studied range, has a negligible impact on the load history maximum values. However, it is important to assess its effect on the interlaminar damage predictions. For this purpose, the output identifier CSQUADSCRT was used to measure the damage initiation in the cohesive surfaces. This variable indicates if the quadratic contact stress damage initiation criterion presented in Equation (9) has been satisfied. When its value reaches 1, damage in the cohesive surface is predicted to initiate. The scalar stiffness degradation for cohesive surfaces, output identifier CSDMG, was used to measure delamination after damage initiation. When it reaches the value of 1, the interface can be considered as fully delaminated (complete debonding). results show that the use of coarser FE meshes ( = 1 mm and = 2 mm) induces hig oscillations on the force-time and force-displacement curves. Nonetheless, the maxim force, maximum displacement, and contact time values are barely affected. This can observed in Table 4, where the values and percentage difference between = 0.3 and the remaining element lengths are presented.   The effect of the FE mesh discretization and mass scaling on the delamination initiation and progression is shown in Figure 10. The results are expressed in terms of the percentage of nodes of the 3D FE mesh for which the CSQUADSCRT is equal to 1 and for those where the CSDMG is higher than 0.6. It can be appreciated that the percentage of nodes where interlaminar damage initiation is predicted decreases for finer FE mesh discretization. If mass scaling is employed, this behavior will be especially obvious. The results indicate different behavior for the fraction of delaminated nodes with and without mass scaling. With mass scaling, the percentage of delaminated nodes slightly increases with the FE mesh refinement but reduces without it. However, for the baseline FE mesh discretization, that is, l e = 0.3 mm, the percentage of delaminated nodes is comparable: 2.8% and 3.22% with mass scaling and without mass scaling, respectively.   Therefore, it is clear that the FE mesh discretization, within the studied rang negligible impact on the load history maximum values. However, it is important t Figure 9. Effect of the FE mesh size on the low-velocity impact response of the force-displacement impact curves. Data suggest that the choice of the FE mesh size does not have a significant impact on the maximum load, maximum displacement, or contact time. However, it affects the interlaminar damage initiation and progression. Consequently, it is recommended to apply the equations suggested by Turon et al. [31] to compute the FE mesh size in order to assure accurate numerical predictions of delamination. The results indicate that the defined semi-automatic mass scaling parameters have no significant impact on the numerical predictions. In addition, taking into consideration the fact that using mass scaling reduces the computational cost of the solutions by around 51%, its application is recommended.
cretization. If mass scaling is employed, this behavior will be especially obvious. The results indicate different behavior for the fraction of delaminated nodes with and without mass scaling. With mass scaling, the percentage of delaminated nodes slightly increases with the FE mesh refinement but reduces without it. However, for the baseline FE mesh discretization, that is, = 0.3 mm, the percentage of delaminated nodes is comparable: 2.8% and 3.22% with mass scaling and without mass scaling, respectively. Data suggest that the choice of the FE mesh size does not have a significant impact on the maximum load, maximum displacement, or contact time. However, it affects the interlaminar damage initiation and progression. Consequently, it is recommended to apply the equations suggested by Turon et al. [31] to compute the FE mesh size in order to assure accurate numerical predictions of delamination. The results indicate that the defined semi-automatic mass scaling parameters have no significant impact on the numerical predictions. In addition, taking into consideration the fact that using mass scaling reduces the computational cost of the solutions by around 51%, its application is recommended.

Numerical-Experimental Correlation
The numerical results obtained with the nine-layer composite laminates FE model with mass scaling and with the different FE mesh sizes are compared with the experimental data presented in [16]. Notice that a higher post-processing by the drop-tower is

Numerical-Experimental Correlation
The numerical results obtained with the nine-layer composite laminates FE model with mass scaling and with the different FE mesh sizes are compared with the experimental data presented in [16]. Notice that a higher post-processing by the drop-tower is responsible for the experimental curves' higher smoothness. The numerical predictions and experimental results are summarized in Table 5 (maximum force, maximum displacement, contact time). It can be observed that the maximum load is slightly overestimated by all the numerical models, presenting errors ranging from 5% to 8.6%, among which the FE model with l e = 0.3 mm was the closest to the experimental averaged value. The maximum displacement is not affected by the FE mesh. Its value is correctly predicted with an error of 2.7% for all FE models. This is due to the fact that the displacement is controlled by the same experimental velocity-time curve that was incorporated to the FE models. The contact time is negligibly affected by the FE mesh size. The numerical predictions are slightly overestimated, presenting an error ranging from 10.4% to 12.3%. Although all FE models present comparable values, the one that best fits the experimental evidence is obtained with l e = 0.3 mm. Moreover, taking also into consideration the results presented in Figure 10, in which it was observed that the FE mesh size considerably affects the interlaminar damage initiation and progression, the use of l e = 0.3 mm, calculated using the equations of Turon et al. [31], is recommended. The numerical and experimental results with l e = 0.3 mm are compared graphically in Figure 11   The numerical and experimental results with = 0.3 mm are compared graphically in Figure 11 (force-time), Figure 12 (energy-time) and Figure 13 (force-displacement). The numerical and experimental curves, during the loading and unloading stages, show a very satisfactory agreement.

Conclusions
Finite element models were generated to study the low-velocity impact response of semicylindrical woven composite laminate shells using ABAQUS/Explicit. These FE models represent a nine-layer laminate with a thickness of 1.6 mm. A continuum damage

Conclusions
Finite element models were generated to study the low-velocity impact response of semicylindrical woven composite laminate shells using ABAQUS/Explicit. These FE models represent a nine-layer laminate with a thickness of 1.6 mm. A continuum damage model at the lamina level to account for intralaminar damage and a surface-based cohe-

Conclusions
Finite element models were generated to study the low-velocity impact response of semicylindrical woven composite laminate shells using ABAQUS/Explicit. These FE models represent a nine-layer laminate with a thickness of 1.6 mm. A continuum damage model at the lamina level to account for intralaminar damage and a surface-based cohesive damage model at the laminas' interface to account for interlaminar damage were used as the two constitutive models to simulate the material damage brought on by low-velocity impact load. The premise for this study was the authors' previous work on the experimental evaluation of the effect of thickness on the impact response of semicylindrical composite laminated shells. As a result, the FE models were tested against the experimental findings for validation purposes. In this way, the stiffness properties of the woven fabric composite laminas were estimated from the constituents' properties of the tested specimens using a homogenization process, while the reaming properties were obtained in the literature. The developed FE models require a fine FE mesh discretization to properly define the interlaminar damage behavior, making them computational expensive. Therefore, numerical simulations were carried out to find an acceptable trade-off between the accuracy of the predictions and the computational cost of the solutions. An efficient approach is employed taking advantage of the model symmetries, continuum shell elements, which have lower computing costs than solid elements, and cohesive surfaces, which do not require element definition. Moreover, the FE mesh discretization and mass-scaling technique were also examined to ascertain how they affect the effectiveness and reliability of the FE model.
The results obtained indicate that the low-velocity impact response of semicylindrical woven composite laminate shells can be reliably predicted using the aforementioned FE models. Furthermore, the mass-scaling technique is successfully used to increase the stable time increment without compromising the accuracy of the dynamic response. Its implementation reduced the computing cost of the solutions by about 51%. The maximum load, maximum displacement, and contact time do not appear to be significantly affected by the choice of FE mesh size. Yet, it has an impact on the initiation and development of interlaminar damage. In order to ensure accurate numerical predictions of delamination, it is advised to use the formula proposed by the literature to compute the characteristic length of the elements of the FE mesh.
The load and energy histories were put to the test in comparison to the experimental findings from low-velocity impacts on specimens with identical elastic properties. A satisfactory correlation between the numerical outcomes and the experimental data was observed. The results show that the maximum load and contact time are slightly overestimated by the FE model, while the maximum displacement is correctly predicted. Nevertheless, all the numerical predictions are comparable with the experimental data.
The numerical simulations complement the previous experimental work, and the developed FE models can be used for further studies, for example, to analyze the response to multi-impacts, the effect of boundary condition or how the geometric parameters of the shells, such as thickness, length, or curvature, affect the impact response of woven fabric composite laminate shells. Acknowledgments: This research was sponsored by national funds through FCT-Fundação para a Ciência e a Tecnologia, under the project UIDB/00285/2020 and LA/P/0112/2020.

Conflicts of Interest:
The authors declare no conflict of interest.